The combined signatures of the tumour microenvironment and nucleotide metabolism-related genes provide a prognostic and therapeutic biomarker for gastric cancer

The tumour microenvironment (TME) is vital to tumour development and influences the immunotherapy response. Abnormal nucleotide metabolism (NM) not only promotes tumour cell proliferation but also inhibits immune responses in the TME. Therefore, this study aimed to determine whether the combined signatures of NM and the TME could better predict the prognosis and treatment response in gastric cancer (GC). 97 NM-related genes and 22 TME cells were evaluated in TCGA-STAD samples, and predictive NM and TME characteristics were determined. Subsequent correlation analysis and single-cell data analysis illustrated a link between NM scores and TME cells. Thereafter, NM and TME characteristics were combined to construct an NM-TME classifier. Patients in the NMlow/TMEhigh group exhibited better clinical outcomes and treatment responses, which could be attributed to the differences in immune cell infiltration, immune checkpoint genes, tumour somatic mutations, immunophenoscore, immunotherapy response rate and proteomap. Additionally, the NMhigh/TMElow group benefited more from Imatinib, Midostaurin and Linsitinib, while patients in the NMlow/TMEhigh group benefited more from Paclitaxel, Methotrexate and Camptothecin. Finally, a highly reliable nomogram was developed. In conclusion, the NM-TME classifier demonstrated a pretreatment predictive value for prognosis and therapeutic responses, which may offer novel strategies for strategizing patients with optimal therapies.


Methods
Data collection. RNA-sequencing (RNA-seq) and the matched clinical characteristics of the TCGA-STAD cohort were obtained from The Cancer Genome Atlas (TCGA) database. Additionally, RNA-seq and clinical data of the GSE84437 cohort, comprising 433 GC samples, were downloaded from the Gene Expression Omnibus (GEO) database as a validation set 24 . Log transformed expression data from raw hybridisation arrays were downloaded and normalised using robust multi-array averaging 25 .
A total of 97 NMRGs were obtained from the Molecular Signatures Database 26 (Supplementary Table S1). For TME cells, we employed CIBERSORT, a deconvolution algorithm 27 , to determine the relative proportions of 22 different types of immune cells. CIBERSORT enrichment values were used to indicate the quantity of each TME cell type in each tumour sample across all cohorts.
Untargeted metabolomic strategies. A total of 33 patients with GC and 27 healthy volunteers were selected and enrolled from the First Affiliated Hospital of Dalian Medical University. The First Affiliated Hospital of Dalian Medical University's institutional ethics committee approved this study. All included patients and healthy volunteers provided consent to the use of their blood samples for research by signing a written consent form. All methods were performed in accordance with the relevant guidelines and regulations. Using gastroscopy or postoperative pathological investigation, the diagnosis of GC was confirmed in patients. Moreover, no evidence of tumours was observed in healthy volunteers. The subjects' clinical features are presented in Supplementary Table S2.
First, 150 μl of each sample was transferred to 1 ml 96-well plates, and 600 μl of methanol was added to precipitate the protein. The mixture was then vortexed for 5 min and centrifuged at 5300 RPM for 20 min (4 °C). Second, two replicates of the 200 μl upper layer were transferred to 450 μl 96-well plates, wherein the samples were concentrated and dried via vacuum centrifugation. These two plates were used for positive and negative ion detection using untargeted metabolomics analysis. Third, the remaining upper layers of all samples were mixed and distributed into 200 μl replicates for use as quality control (QC) samples. Finally, polar metabolite analysis was performed on an Ultimate 3000 ultra-high-performance liquid chromatograph and Q-Orbitrap mass spectrometer.
Establishment of NM score, TME score and NM-TME classifier. The differentially expressed NMRGs between GC and normal tissues were identified using the 'limma' packages (P < 0.05) 28 . Following this, differentially expressed genes (DEGs) were evaluated using univariate Cox regression analysis to acquire the genes with prognostic value (P < 0.05). Meanwhile, genes with prognostic significance were validated using Kaplan-Meier (KM) analysis 29 . These genes were then incorporated into a multivariate Cox regression model to calculate NM scores using the following formula: NM scores = n i=1 Expi * Coefi (where n, Coefi and Expi denote the number of prognostic genes, the expression value and the coefficient of gene i, respectively). As for TME cells, 22 immune cell enrichment scores in patients with GC were calculated and then immune cells with prognostic significance were identified using KM analysis. Additionally, the TME score was calculated using the following formula: n i=1 Expi * Coefi (where n, Coefi and Expi denote the number of prognostic immune cell, the infiltrating value and the coefficient of immune cell i, respectively). Gene set expression analysis (GSEA) was used to analyze the potential functions of different groups in TCGA data set. Signaling pathway differences were integrated through Kyoto Encyclopedia of Genes and Genomes (KEGG) detabase [30][31][32] . We then integrated NM and TME scores to construct an NM-TME classifier and classified tumours into the following subgroups: NMlow/TMEhigh, Single-cell data processing and cell-cell communication. We next analysed the relationship between NM scores and immune cells using single-cell data. GSE167297 was used to obtain single-cell transcriptome data from 10 GC samples 34 . Seurat R package was used to analyse the single-cell RNA-seq data 35 . Additionally, the 'NormalizeData' and 'FindVariableFeatures' functions in the Seurat package were used to normalise the count and expel cells containing less than 200 genes, more than 2500 genes, more than 20% of mitochondria or more than 3% of haemoglobin and then identify the 3000 highly variable genes. Moreover, the non-linear dimensional reduction was performed using the UMAP and tSNE methods. Cluster biomarkers were identified using the 'Seurat' package. The 'CellChat' R package's method of identifying communication molecules at single-cell resolution was used to analyse the relationships between cells that were involved in communication 36 .
Weighted gene co-expression network analysis (WGCNA). To further investigate the potential reasons for the significant differences between the NMlow/TMEhigh and NMhigh/TMElow groups, we performed WGCNA using the R package 'WGCNA' . A weighted value was set that conformed to the scale-free network law (scale-free R 2 = 0.9). Topological coefficients were employed to determine the degree of dissimilarity between nodes and create a hierarchical clustering tree to separate modules. The modules with the highest correlation to the NMlow/TMEhigh and NMhigh/TMElow groups were considered key modules. Finally, functional enrichment analysis of key module genes was performed using the Metascape database 37 .

Immune cell infiltration and immune checkpoint gene (ICG) expression between different groups.
Tracking tumour immunophenotype (TIP) is a database that aids in the understanding of the mechanism of tumour immune activity and the proportion of immune cell infiltration 38 . TIP follows a seven-step 'cancer-immunity cycle' analysis 39 , wherein stepwise events are grouped into 23 groups with 178 stimulatory or inhibitory signature genes. Herein, TIP was used to explore the immune cell infiltration between different groups. We also examined the differential expression of common ICGs between the different groups, showing only statistically significant results.
Tumour mutation burden (TMB) and immunotherapy response analysis. In this study, mutation data were downloaded from the TCGA database. The top 20 most frequently mutated genes in different groups were identified using the 'maftools' package in R. Based on the somatic mutation data in each tumour, TMB was calculated as the number of mutated bases per million bases and compared across groups. Subsequently, we also explored the survival probability between different TMB and risk scores to highlight the crucial role of TMB in GC. To further estimate the response of immunotherapy, the tumour immune dysfunction and exclusion (TIDE) algorithm was used 40 . Additionally, the Immune Checkpoint Inhibitor (ICI) Immunophenoscore (IPS) file from The Cancer Immunome Atlas Database was retrived 41 . The immunotherapeutic relevance of the signature was evaluated using IPS, a reliable tool for assessing tumour immunogenicity.

Nomogram construction and validation.
To determine if the risk score was an independent predictor of GC, univariate and multivariate Cox regression analyses were performed. Clinicopathological parameters and risk scores were considered in the development of the nomogram model for predicting the prognosis of GC using the 'rms' R package 42 . Additionally, we used a ROC curve to assess the validity of the established nomogram.
Statistical analyses. The statistical analyses were conducted using strawberry-Perl and R software (R-4.13).

Results
Landscape of the genetic variation of NMRGs in GC. Figure 1 presents the workflow of the study.
The prognostic values of NM and TME score. To construct an NM prognostic model, we first performed a univariate Cox survival analysis on 77 differentially expressed NMRGs, of which six were statistically significant (Supplementary Table S3). Additionally, the prognostic significance of the six genes was validated using KM analysis ( Supplementary Fig. 3A). Furthermore, a heatmap of the expression of the six genes in tumour and normal tissues was also drawn ( Fig. 2A). We then subjected the six genes to multivariate cox analysis (Fig. 2B) and correlation coefficients were calculated to construct a model (Supplementary Table S4). The NM score was calculated for each patient, and the patients were classified into high and low score groups based on the median value. The KM curve showed that the high-risk patients had a worse prognosis (Fig. 2C).
Regarding the TME prognostic model, a high infiltration of activated CD4 memory-activated T cells, CD8 T cells and activated dendritic cells (DCs) were observed to be associated with a better prognosis for patients with GC ( Supplementary Fig. 3B). Similarly, these cells were subjected to multivariate cox analysis (Fig. 2D) and correlation coefficients were calculated to construct a model (Supplementary Table S5). The KM curve showed that high-TME score samples had a better survival prognosis than those with low-TME scores (Fig. 2E). GSEA revealed that the high NM score group was mainly enriched in cancer-related and classical oncogenic pathways, while the high TME score group was mainly enriched in immune-related pathways. (Supplementary Fig. 3C,D).
Single-cell data analysis to explore the association between NM scores and TME cells. First, we investigated the correlation between the six NM model genes and the three TME cells. We found that T cells www.nature.com/scientificreports/ CD8 were negatively correlated with UPP1, ENTPD2, NT5E and positively correlated with DPYS and AK1; T cells CD4 memory activated were negatively correlated with AK5, ENTPD2, NT5E, DPYS, AK1; dendritic cells were negatively correlated with AK5, ENTPD2, DPYS, and positively correlated with UPP1 (Fig. 3A). To further explore their association, we downloaded single-cell data from the GEO database, comprising 10 GC samples. The clustering and annotated results are presented in Fig. 3B. Subsequently, we calculated the NM scores in different cell types and found that the NM scores were significantly higher in monocytes and endothelial cells than in B cells, T cells, CD8+ T cells, epithelial, macrophages, Tregs and mast cells (Fig. 3C,D). Based on the NM score, monocytes and endothelial cells were divided into low NM score, medium NM score and high NM score monocytes and endothelial cells for cell communication analysis. The monocytes and endothelial cells with low NM scores had more abundant communication with other immune cells (Fig. 3E-H). Therefore, low NM scores could have a synergistic effect with high TME scores and combining the NM model with the TME model may be a feasible method.
NM-TME classifier construction and validation. Next, we constructed the NM-TME classifier by combining the NM and TME scores. It divided patients with GC into four categories: NMhigh/TMEhigh, NMhigh/ TMElow, NMlow/TMEhigh and NMlow/TMElow. Survival analysis revealed that the NMhigh/TMElow group had a poorer prognosis while the NMlow/TMEhigh group had a better prognosis among the groups (Fig. 4A). Patients in the NMhigh/TME high and NMlow/TME low subgroups showed less divergent prognoses. As a result, we combined them to form a mixed subgroup (Fig. 4B). Additionally, the area under the curve (AUC) values of the NM-TME classifier were 0.732, 0.708, 0.702 and 0.807 for 1, 3, 5 and 7 years, respectively (Fig. 4C), indicating that the NM-TME classifier plays a significant role in the survival prediction of patients with GC. Furthermore, we also verified the prognostic significance of the NM-TME classifier in the GEO cohort, which revealed significant prognostic differences between the groups (Supplementary Fig. 4A). Moreover, the  Functional enrichment analysis and WGCNA. Functional enrichment of the three groups revealed that the NMhigh/TMElow group was mainly enriched in the regulation of the olefinic compound metabolic process, endothelial cell differentiation and stem cell proliferation, while the NMlow/TMEhigh was majorly positively associated with the positive regulation of T cell migration and negatively associated with the canonical Wnt signalling pathway (Fig. 4D). Furthermore, WGCNA identified four modules (Fig. 5A,B). Among them, the turquoise module was most relevant and opposite to each other for the NMlow/TMEhigh and NMhigh/TMElow groups. Therefore, the turquoise module gene could be associated with significantly different prognoses between the NMlow/TMEhigh and NMhigh/TMElow groups. Using the Metascape database, enrichment analysis of these genes revealed that they were mainly enriched in vasculature development, NABA core matrisome and extracellular matrix organization (Fig. 5C).

Differences in immune cell infiltration and ICG expression based on the NM-TME classifier.
First, we compared the abundance of immune cell infiltration between the different groups. The immune www.nature.com/scientificreports/ cell infiltration was more abundant in the NMlow/TMEhigh group, especially CD8 T cells, Th1 cells, NK cells, CD4 T cells and macrophages (Fig. 6A). Notably, the better prognosis in the NMlow/TMEhigh group could be attributed to the abundant immune cell infiltration. Meanwhile, we also explored whether the expression of common ICGs differed between the groups. Most ICGs were differentially expressed between the groups, with high expression observed in the NMlow/TMEhigh group (Fig. 6B). These differentially expressed ICGs could be potential therapeutic targets. Additionally, it also suggests that NMlow/TMEhigh patients may benefit more from immune checkpoint blockade (ICB) therapy. HLA is a polygenic and polymorphic complex involved in antigen presentation 43 . Figure 6C shows that HLA-B, HLA-C, HLA-F and HLA-DOB were expressed the highest in the NMlow/TMEhigh group.

Intergroup differences in cancer somatic mutations.
Numerous studies have demonstrated the association between somatic mutations in tumour genomes and the response to immunotherapy 44 . We therefore examined the TMB distributions among the various groups based on the NM-TME classifier. The NMlow/TMEhigh group had a higher TMB, while the NMhigh/TMElow group had a lower TMB, indicating that the NMlow/ TMEhigh group may benefit more from immunotherapy (Fig. 7A). Additionally, the NMhigh/TMElow/TMBhigh group had a lower prognosis than patients in the other groups (Fig. 7B). Figure 7C,D display the top 20 genes with high mutation frequencies in the NMlow/TMEhigh and NMhigh/TMElow groups.
Personalised treatment based on the NM-TME classifier. Considering that drugs targeting PD-1 and CTLA-4 have recently received approval for the treatment of several cancers, we evaluated whether the NM-TME classifier could predict patients' reactions to immunotherapy. The patients in the NMlow/TMEhigh group were observed to have a better response rate to immunotherapy than the other two groups (Fig. 8A). Microsatellite instability-high (MSI-H) is a potential predictor of immunotherapy response targeting PD-1 or its ligand PD-L1 45 . Accordingly, the proportion of MSI-H in the NMlow/TMEhigh group was higher than that in the other two groups (Fig. 8B). Additionally, we investigated the relationship between the NM-TME classifier and IPS in patients with GC to predict the response to ICIs. Figure 8C-F presents the differences in the results of CTLA-4/ www.nature.com/scientificreports/ PD-1 inhibitor treatment between the NMlow/TMEhigh and Nmhigh/TMElow groups. The NMlow/TMEhigh group has higher IPS scores, implying more immunogenicity in the NMlow/TMEhigh group. Furthermore, we performed a difference analysis between the immunotherapy-responsive and non-responsive groups and also the NMlow/TMEhigh and NMhigh/TMElow groups. DEGs were then analysed using the Proteomaps 2.0 database 46 . Notably, the pattern of proteomap in the NMlow/TMEhigh group and immunotherapy-responsive groups were similar (Fig. 8G,H). These findings suggest that the NM-TME classifier can be used to predict patients' responses to immunotherapy. Given that targeted therapy is an effective approach in the treatment of GC, it has important clinical applications and prospects. We, therefore, investigated whether the NM-TME classifier could predict drug sensitivity in patients with GC. The NMhigh/TMElow group benefited more from Imatinib, Midostaurin and OSI-906 (Linsitinib), while those in the NMlow/TMEhigh group benefited more from Paclitaxel, Methotrexate and Camptothecin ( Supplementary Fig. 5A-F).

Nomogram development and verification. Univariate and multivariate Cox regression analyses indi-
cated that the NM-TME classifier was an independent predictor of prognosis with the highest hazard ratio (HR) (Fig. 9A,B). Following this, the NM-TME classifier and clinical features were combined to construct a nomogram. To predict the survival of patients with GC over 1 to 5 years, the values of each variable can be added to obtain the total score (Fig. 9C). Moreover, the AUC values of the nomogram for 1-, 3-and 5-year OS were 0.826, 0.841 and 0.822, respectively (Fig. 9D).

Discussion
Owing to its high morbidity, the poor incidence of early diagnosis and low survival rate, GC poses a severe threat to the populations worldwide 47 . Moreover, the development of tumours is consistent with abnormal metabolism 48 . Recent studies have demonstrated that aberrant NM speeds up the progression of tumours while suppressing the TME's normal immune response 49,50 . Therefore, to treat malignancies and prevent recurrence and metastasis, the intervention or regulation of molecular pathways related to aberrant NM in malignant cells has emerged as a novel therapeutic strategy 20 . The TME has a vital role in tumour development, growth, metastasis and therapeutic response 51 . As a therapeutic target in tumours, TME has attracted significant research and clinical interest 52 . However, very few studies have reported on the use of NM-TME characteristics in predicting GC prognosis and treatment response. In this study, we combined the NM and TME features, for the first time, to construct the NM-TME classifier, which consists of three types of immune cells and six NMRGs. Currently, in various malignancies, next-generation sequencing is becoming a complementary diagnostic tool that guides decision-making to achieve precise and personalised therapy regimens. Accordingly, based on our constructed www.nature.com/scientificreports/ NM-TME signature, clinicians can quantify specific immune cells and target genes using sequencing technology and related algorithms to effectively predict prognosis, immunotherapy response and targeted therapy response in patients with GC. First, we constructed an NM prognostic model with six genes, namely AK1, DPYS, NT5E, ENTPD2, AK5 and UPP1. Second, we constructed a prognostic model of the TME, consisting of activated DCs, activated CD4 memory T cells and CD8 T cells. Both prognostic models classified patients with GC into two groups with significant prognostic differences. Additionally, patients with high NM scores had a poor prognosis and majorly played a role in cancer-related pathways. In contrast, patients with high TME scores had a better prognosis and were mainly involved in immune-related pathways. Therefore, we speculated that these two models could have synergistic effects.
Moreover, a correlation was observed between the six NMRGs and three TME cells that were involved in the construction of the model. Additionally, we used single-cell data to further explore the association between NM scores and immune cells. After QC, clustering, and annotation of the single cell data, we calculated the NM scores in each cell type. The results showed that the NM scores in monocytes and endothelial cells were significantly higher than in other cells. We then divided monocytes and endothelial cells into monocytes and endothelial cells with high NM, medium NM and low NM scores. Furthermore, cell communication analysis also showed that low NM score monocytes and endothelial cells were more closely related to other immune cell types. Thus, these findings suggest a strong association between low NM scores and high TME scores, highlighting their synergistic effect on the prognosis of patients with GC.
Based on the above analysis, we constructed an NM-TME classifier that can classify patients with GC into different subgroups based on NM and TME scores. Survival analysis showed significant differences in prognosis between the groups, which was consistent with the results of the test set. A key module was identified to be significantly associated with the NMlow/TMEhigh and NMhigh/TMElow groups via WGCNA, which could be responsible for their significant differences. Moreover, the key module genes were mainly enriched in vasculature development, NABA core matrisome and extracellular matrix organization.
We then examined the immune status of the different groups of patients based on the NM-TME classifier. Patients in the NMlow/TMEhigh group had a higher abundance of immune cell infiltration, which may be the reason for better prognosis in NMlow/TMEhigh group. Moreover, ICB therapy as emerging immunotherapy target has demonstrated therapeutic efficacy in the treatment of human malignancies 53 . Herein, most ICGs were www.nature.com/scientificreports/ highly expressed in the NMlow/TMEhigh group. These differentially expressed ICGs could be potential therapeutic targets, suggesting that patients in the NMlow/TMEhigh group may benefit more from ICB. TMB has been demonstrated to be utilised as a predictor of ICB efficacy and has become a biomarker in certain types of cancer to identify patients who might benefit from immunotherapy 44,54,55 . On analysing TMB values, the NMlow/TMEhigh group exhibited a higher TMB, while a converse trend was observed in the NMhigh/ TMElow group. Thus, patients in the NMlow/TMEhigh group can be considered more sensitive to ICB treatment.
According to recent studies, blocking PD-1 is not inferior to chemotherapy 56 and combining a PD-1 inhibitor with chemotherapy improves survival in individuals with advanced GC compared to chemotherapy alone 57 . Nevertheless, anti-PD-1 immunotherapy has been reported to be efficient in 15-60% of patients. Therefore, we investigated the relationship between the NM-TME classifier and the outcome of CTLA-4/PD-1 inhibitor therapy. The NMhigh/TMElow group had a better response to CTLA-4/PD-1 inhibitor therapy. Furthermore, validation using the TIDE database revealed that the immunotherapy response rate and proportion of MSI-H were higher in the NMlow/TMEhigh group. Meanwhile, the pattern of proteomap in the NMlow/TMEhigh group and immunotherapy responder were also similar. These results further demonstrate that patients in the NMlow/TMEhigh group are more sensitive to immunotherapy and that the NM-TME classifier can effectively predict patients' response to immunotherapy.
In addition, the NM-TME classifier was able to predict chemotherapy drug sensitivity. The NMhigh/TMElow group benefited more from Imatinib, Midostaurin and OSI-906 (Linsitinib), while patients in the NMlow/TMEhigh group benefited more from Paclitaxel, Methotrexate and Camptothecin. Drug resistance is a major challenge in cancer treatment. It is speculated that drug resistance in cancer is driven by genetic mutations. Despite the unclear mechanism of drug resistance, there is evidence for an important role of reversible proteomic and epigenetic mechanisms in drug resistance. Additionally, mechanisms mediated by the TME and tumour heterogeneity greatly contribute to cancer therapy resistance 58 . Tyrosine kinase inhibitors (TKIs) therapy, such as Imatinib and OSI-906 (Linsitinib), play a role in the TME remodelling and enhance therapeutic response, but TME changes can also induce drug resistance and promote tumour growth. The higher resistance to Imatinib and OSI-906 (Linsitinib) in patients in the NMlow/TMEhigh group could be attributed to their more abundant immunosuppressive TMEs, such as macrophages, neutrophils, Tregs and myeloid-derived suppressor cells (MDSCs) 59 . It has been reported that Midostaurin may enhance anti-tumour effects by modulating the distribution of immune cells in the TME. Thus, resistance to Midostaurin could also be associated with the anti-tumour immunity of neutrophils and MDSCs 60 . Paclitaxel can promote the polarisation to DCs and the proliferation and activity of CD8+ T cells and NK cells to exert stronger anti-tumour effects, which leads to a higher sensitivity to paclitaxel www.nature.com/scientificreports/ in the NMlow/TMEhigh group 61 . Methotrexate is an anti-tumour agent that interferes with folic acid metabolism. Studies have reported that higher NADPH levels in acute myeloid leukaemia promote Methotrexate resistance and that NADPH is involved in nucleotide synthesis 62 , which could be associated with greater sensitivity in the NMlow/TMEhigh group. However, the corresponding mechanisms involved in GC require further investigation. The multifactorial analysis demonstrated that the NM-TME classifier is an independent prognostic factor for patients with GC, with excellent prognostic predictive power. Finally, to fully exploit the prognostic potential of the NM-TME classifier, the survival rate of patients with GC was quantified after constructing a nomogram based on the signature and clinical features. The ROC curve illustrated the high-precision predictive capability of the nomogram.
To our knowledge, this is the first study to use bioinformatics to combine TME and NM features to analyse their role in GC prognosis, immunotherapy and chemotherapy. Nonetheless, this investigation is not without its drawbacks. The data used herein are from online databases, namely TCGA and GEO. Thus, these findings require further validation using real prospective clinical cohorts. Furthermore, basic investigations on the function of the TME and NM in the aetiology and progression of GC are required as the current understanding of this topic is limited.

Conclusion
In our study, an NM-TME signature was constructed by combining NM and TME features to predict the prognosis, immunotherapy and chemotherapy effects of patients with GC. This classifier has been well-validated from different points of view.

Data availability
The datasets analysed in this study are publicly available from the TCGA and GEO (GSE84437, https:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE84 437) databases. Furthermore, the raw data and analytic technologies used in this study can be obtained from the corresponding author and first author upon reasonable request.